/*
 * Copyright (c) 2009, Luke Emrose
 *
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * Redistributions of source code must retain the above copyright
 * notice, this list of conditions and the following disclaimer.
 * Redistributions in binary form must reproduce the above copyright
 * notice, this list of conditions and the following disclaimer in the
 * documentation and/or other materials provided with the
 * distribution.  Neither the name of Luke Emrose Sole Trader nor the
 * names of its contributors may be used to endorse or promote
 * products derived from this software without specific prior written
 * permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE
 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
 * OF THE POSSIBILITY OF SUCH DAMAGE.
 */
 
 #ifndef __WAVELET_NOISE_H__
#define __WAVELET_NOISE_H__

/**
 * @file
 * Reference: Wavelet Noise - Robert L. Cook, Tony DeRose, Pixar Animation Studios
 *
 * Original Note:
 * this code is designed for brevity, not efficiency; many operations can be hoisted,
 * precomputed, or vectorized. Some of the straightforward details, such as tile meshing,
 * decorrelating bands and fading out the last band, are omitted in the interest of space.
 *
 * Modified Code:
 * Luke Emrose - evolutionarytheory@gmail.com
 */

#include "noisesGlobals.hpp"

// std includes
#include <cmath>

// boost includes
#include <boost/random.hpp>

// numerical precision to calculate noise with
typedef RtFloat Precision;

static Precision *noiseTileData;
static int noiseTileSize;

#define ARAD 16

// define the type of random number generator to use
typedef boost::minstd_rand base_generator_type;

// define the random number generator state variables
base_generator_type generator( base_generator_type( 42u ) );
boost::normal_distribution< Precision > norm_dist( boost::normal_distribution< Precision >( 0, 1 ) );
boost::variate_generator< base_generator_type&, boost::normal_distribution< Precision > > norm( boost::variate_generator< base_generator_type&, boost::normal_distribution< Precision > >( generator, norm_dist ) );

/*
 * @brief
 * Modulus function
 */
template < typename T >
T Mod( T x, T n )
{
	T m = x % n;
	return ( m < 0 ) ? m + n : m;
}

/*
 * @brief
 * Downsample a noise tile
 */
template < typename T >
void Downsample
(
	T *from,
	T *to,
	int n, int
	stride
)
{
	T *a, aCoeffs[ 2 * ARAD ] =
	{
		0.000334f,	-0.001528f,	0.000410f,	0.003545f,	-0.000938f,	-0.008233f,	0.002172f,	0.019120f,
		-0.005040f,	-0.044412f,	0.011655f,	0.103311f,	-0.025936f,	-0.243780f,	0.033979f,	0.655340f,
		0.655340f,	0.033979f,	-0.243780f,	-0.025936f,	0.103311f,	0.011655f,	-0.044412f,	-0.005040f,
		0.019120f,	0.002172f,	-0.008233f,	-0.000938f,	0.003546f,	0.000410f,	-0.001528f,	0.000334f
	};

	a = &aCoeffs[ARAD];

	for( int i = 0; i < n / 2; ++i )
	{
		to[ i * stride ] = 0;
		for( int k = 2 * i - ARAD; k <= 2 * i + ARAD; ++k )
		{
			to[ i * stride ] += a[ k - 2 * i ] * from[ Mod( k, n ) * stride ];
		}
	}
}

/*
 * @brief
 * Upsample a noise tile
 */
template < typename T >
void Upsample
(
	T *from,
	T *to,
	int n,
	int stride
)
{
	T *p, pCoeffs[ 4 ] = { 0.25, 0.75, 0.75, 0.25 };

	p = &pCoeffs[2];

	for( int i = 0; i < n; ++i )
	{
		to[ i * stride ] = 0;
		for( int k = i / 2; k <= i / 2 + 1; ++k )
		{
			to[ i * stride] += p[ i - 2 * k ] * from[ Mod( k, n / 2 ) * stride ];
		}
	}
}

/*
 * @brief
 * Generate a noise tile
 */
template < typename T >
void initWaveletNoiseTile
(
	int n,
	int olap
)
{
	if( n % 2 )
	{
		n++; /* tile size must be even */
	}

	int ix, iy, iz, i;
	int sz = n * n * n;

	T *temp1 = new T[ sz ];
	T *temp2 = new T[ sz ];
	noiseTileData = new T[ sz ];

	// Step 1. Fill the tile with random numbers in the range -1 to 1.
	for( i = 0; i < n * n * n; ++i )
	{
		// insert guassian noise in the range -1 to 1 here....
		//gaussianNoise< float >();
		noiseTileData[ i ] = norm() * ( T ) 0.5 + ( T ) 0.5;
		//noiseTileData[ i ] = ( T ) rand() / ( T ) RAND_MAX;
	}

	// Steps 2 and 3. Downsample and upsample the tile
	for( iy = 0; iy < n; ++iy )
	{ 
		for( iz = 0; iz < n; ++iz )
		{
			// each x row
			i = iy * n + iz * n * n;
			Downsample( &noiseTileData[ i ], &temp1[ i ], n, 1 );
			Upsample( &temp1[ i ], &temp2[ i ], n, 1 );
		}
	}

	for( ix = 0; ix < n; ++ix )
	{
		for( iz = 0; iz < n; ++iz )
		{
			// each y row
			i = ix + iz * n * n;
			Downsample( &temp2[ i ], &temp1[ i ], n, n );
			Upsample( &temp1[ i ], &temp2[ i ], n, n );
		}
	}

	for( ix = 0; ix < n; ++ix )
	{
		for( iy = 0; iy < n; ++iy )
		{
			// each z row
			i = ix + iy * n;
			Downsample( &temp2[ i ], &temp1[ i ], n, n * n );
			Upsample( &temp1[ i ], &temp2[ i ], n, n * n );
		}
	}

	// Step 4. Subtract out the coarse-scale contribution
	for( i = 0; i < n * n * n; ++i )
	{
		noiseTileData[ i ] -= temp2[ i ];
	}

	// Avoid even/odd variance difference by adding odd-offset version of noise to itself.
	int offset = n / 2;

	if( offset % 2 == 0 )
	{
		offset++;
	}

	for( i = 0, ix = 0; ix < n; ++ix )
	{
		for( iy = 0; iy < n; ++iy )
		{
			for( iz = 0; iz < n; ++iz )
			{
				temp1[ i++ ] = noiseTileData[ Mod( ix + offset, n ) + Mod( iy + offset, n ) * n + Mod( iz + offset, n ) * n * n ];
			}
		}
	}

	for( i = 0; i < n * n * n; ++i )
	{
		noiseTileData[ i ] += temp1[ i ];
	}

	noiseTileSize = n;

	delete [] temp1;
	delete [] temp2;
}

/**
 * @brief
 * clear memory used by wavelet noise tile
 */
void clearWaveletNoiseTile( void )
{
	delete [] noiseTileData;
}

/*
 * @brief
 * 3D wavelet noise function
 */
template < typename T >
T waveletNoise( T p[ 3 ] )
{
	// Non-projected 3D noise
	int i, f[ 3 ], c[ 3 ], mid[ 3 ], n = noiseTileSize; // f, c = filter, noise coeff indices
	T w[ 3 ][ 3 ], t, result = 0;

	// Evaluate quadratic B-spline basis functions
	for( i = 0; i < 3; ++i )
	{
		mid[ i ] = ( int ) ceil( p[ i ] - ( T ) 0.5 );
		t = mid[ i] - ( p[ i ] - ( T ) 0.5 );
		w[ i ][ 0 ] = t * t / 2;
		w[ i ][ 2 ] = ( 1 - t ) * ( 1 - t ) / 2;
		w[ i ][ 1 ] = 1 - w[ i ][ 0 ] - w[ i ][ 2 ];
	}

	// Evaluate noise by weighting noise coefficients by basis function values
	for( f[ 2 ] =- 1; f[ 2 ] <= 1; ++f[ 2 ] )
	{
		for( f[ 1 ] =- 1; f[ 1 ] <= 1; ++f[ 1 ] )
		{
			for( f[ 0 ] =- 1; f[ 0 ] <= 1; ++f[ 0 ] )
			{
				T weight = 1;
				for( i = 0; i < 3; i++ )
				{
					c[ i ] = Mod( mid[ i ] + f[ i ], n );
					weight *= w[ i ][ f[ i ] + 1 ];
				}
				result += weight * noiseTileData[ c[ 2 ] * n * n + c[ 1 ] * n + c[ 0 ] ];
			}
		}
	}
	return result;
}

/*
 * @brief
 * Projected 3D wavelet noise
 */
template < typename T >
T waveletProjectedNoise( T p[ 3 ], T normal[ 3 ] )
{
	// 3D noise projected onto 2D
	int i, c[ 3 ], min[ 3 ], max[ 3 ], n = noiseTileSize; // c = noise coeff location
	T support, result = 0;

	// Bound the support of the basis functions for this projection direction
	for( i = 0; i < 3; ++i )
	{
		support = 3 * abs( normal[ i ] ) + 3 * sqrt( ( 1 - normal[ i ] * normal[ i ] ) / 2 );
		min[ i ] = ( int ) ceil( p[ i ] - ( 3 * abs( normal[ i ] ) + 3 * sqrt( ( 1 - normal[ i ] * normal[ i ]) / 2 ) ) );
		max[ i ] = ( int ) floor( p[ i ] + ( 3 * abs( normal[ i ] ) + 3 * sqrt( ( 1 - normal[ i ] * normal[ i ] ) / 2 ) ) );
	}

	/* Loop over the noise coefficients within the bound. */
	for( c[ 2 ] = min[ 2 ]; c[ 2 ] <= max[ 2 ]; ++c[ 2 ] )
	{
		for( c[ 1 ] = min[ 1 ]; c[ 1 ] <= max[ 1 ]; ++c[ 1 ] )
		{
			for( c[ 0 ] = min[ 0 ]; c[ 0 ] <= max[ 0 ]; ++c[ 0 ] )
			{
				T t, t1, t2, t3, dot = 0, weight = 1;
				// Dot the normal with the vector from c to p
				for( i = 0; i < 3; ++i )
				{
					dot += normal[ i ] * ( p[ i ] - c[ i ] );
				}

				/* Evaluate the basis function at c moved halfway to p along the normal. */
				for( i = 0; i < 3; i++ )
				{
					t = ( c[ i ] + normal[ i ] * dot / 2 ) - ( p[ i ] - 1.5f );
					t1 = t - 1;
					t2 = 2 - t;
					t3 = 3 - t;
					weight *= ( t <= 0 || t >= 3 ) ? 0 : ( t < 1 ) ? t * t / 2 : ( t < 2 ) ? 1 - ( t1 * t1 + t2 * t2 ) / 2 : t3 * t3 / 2;
				}

				/* Evaluate noise by weighting noise coefficients by basis function values. */
				result += weight * noiseTileData[ Mod( c[ 2 ], n ) * n * n + Mod( c[ 1 ], n ) * n + Mod( c[ 0 ], n ) ];
			}
		}
	}

	return result;
}

/*
 * @brief
 * Multiband wavelet noise
 */
template < typename T >
T waveletMultibandNoise( T p[ 3 ], T s, T *normal, int firstBand, int nbands, T *w )
{
	T q[ 3 ];
	T result = 0;
	T variance = 0;
	int i, b;

	for( b = 0; b < nbands && s + firstBand + b < 0; ++b )
	{
		for( i = 0; i <= 2; ++i )
		{
			q[ i ] = 2 * p[ i ] * powf( 2, ( T ) firstBand + b );
		}

		result += ( normal ) ? w[ b ] * waveletProjectedNoise< T >( q, normal ) : w[ b ] * waveletNoise< T >( q );
	}

	for( b = 0; b < nbands; ++b )
	{
		variance += w[ b ] * w[ b ];
	}

	// Adjust the noise so it has a variance of 1.
	if( variance )
	{
		result /= sqrt( variance * ( ( normal ) ? ( T ) 0.296 : ( T ) 0.210 ) );
	}

	return result;
}

#endif